4. Wave Equation#
We consider the wave equation as typical example of second order hyperbolic equation. We have a second order derivative in time, and second order derivatives in space, but of different signs. This completely changes the behaviour of the solution.
Since it is second order in time, we need two initial conditions, for the function itself and for the time derivative:
We use a variational formulation in space and obtain the second order ordinary differential equation (ODE) in time:
with initial conditions
4.1. Newmark time-stepping method#
We reduce the second order ODE to a first order system. For this, we introduce the velocity function \(v = \dot u\):
Now we apply the trapezoidal rule for time-discretization:
Expressing \(u^j\) from the first equation, and inserting it into the second we obtain
and shifting the unknowns to the left we obtain a linear system for the update of velocity:
Although by introducing the system of ODEs, the size of the linear system to solve did not increase.
from ngsolve import *
from ngsolve.webgui import Draw
from time import sleep
from netgen.occ import *
rect = MoveTo(-1,-1).Rectangle (2,2).Face()
rect.edges.Min(X).name="left"
circ = Circle ( (0.7, 0.0), 0.1).Face()
shape = rect-circ
mesh = Mesh (OCCGeometry( shape, dim=2).GenerateMesh(maxh=0.03)).Curve(3)
Draw (shape)
Draw (mesh);
tau = 0.002
tend = 4
u0 = exp(-400*( (x)**2 + (y)**2))
v0 = 0
fes = H1(mesh, order=3)
u,v = fes.TnT()
mform = u*v*dx
aform = grad(u)*grad(v)*dx
m = BilinearForm(mform).Assemble()
a = BilinearForm(aform).Assemble()
mstar = BilinearForm(mform+tau**2/4*aform).Assemble()
mstarinv = mstar.mat.Inverse(inverse="sparsecholesky")
f = LinearForm(fes).Assemble()
gfu = GridFunction(fes)
gfv = GridFunction(fes)
gfu.Set(u0)
gfv.Set(u0)
scene = Draw (gfu, deformation=True, euler_angles=(-60,0,-10))
sleep (5)
with TaskManager():
for j in range(int(tend/tau)):
gfu.vec.data += tau/2 * gfv.vec
gfv.vec.data += tau * mstarinv * (f.vec - a.mat * gfu.vec)
gfu.vec.data += tau/2 * gfv.vec
if j%10 == 0:
scene.Redraw()
scene.Redraw()